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Abstract: We present a detailed description of the methods used 
in our model of mode instability in high-power, rare earth-doped, 
large-mode-area fiber amplifiers. Our model assumes steady-periodic 
behavior, so it is appropriate to operation after turn on transients 
have dissipated. It can be adapted to transient cases as well. We 
describe our algorithm, which includes propagation of the signal 
field by fast-Fourier transforms, steady-state solutions of the laser 
gain equations, and two methods of solving the time-dependent heat 
equation: altcrnating-direction-implicit integration, and the Green's 
function method for steady-periodic heating. 

© 2013 Optical Society of America 

OCIS codes: (060.2320) Fiber optics amplifiers and oscillators; (060.4370) Non- 
linear optics, fibers; (140.6810) Thermal effects; (190.2640) Stimulated scattering, 
modulation, etc. 



References and links 

1. A. V. Smith and J. J. Smith, "Mode instability in high power fiber amplifiers," Opt. Express 19, 
10180-10192 (2011), [http : //www. opticsinf obase . org/oe/abstract . cf m?URI=oe- 19-11-10180 1 

2. A. V. Smith and J. J. Smith, "Influence of pump and seed modulation on the 
mode instability thresholds of fiber amplifiers," Opt. Express 20, 24545-24558 (2012), 
http : //www. opticsinf obase . org/oe/abstract . cf m?URI=oe-20-22-24545 

3. H.-J. Otto, F. Stutzki, F. Jansen, T. Eidam, C. Jauregui, J. Limpcrt, and A. Tiinnermann, 
"Temporal dynamics of mode instabilities in high-power fiber lasers and amplifiers," Opt. Express 
20, 15710-15722 (2012). [http : //www ■ opticsexpress ■ org/abstract . cf m?URI=oe-20- 14- 15710 1 

4. K. R. Hansen. T. T. Alkcskjold, J. Brocng, and J. Laegsgaard, "Thermally induced 
mode coupling in rare-earth doped fiber amplifiers," Opt. Lett. 37, 2382—2384 (2012), 
http : //www. opticsinf obase . org/oe/ abstract . cf m?URI=ol-37- 12-2382 

5. B. Ward, C. Robin, and I. Dajani, "Origin of thermal modal instabilities 
in large mode area fiber amplifiers," Opt. Express 20, 11407-11422 (2012), 
http : //www. opticsinf obase . org/oe/ abstract . cf m?URI=oe-20- 10- 11407 

6. C. Jauregui, T. Eidam, H.-J. Otto, F. Stutzki, F. Jansen, J. Limpcrt, and A. Tiinnermann, 
"Physical origin of mode instabilities in high-power fiber laser systems," Opt. Express 20, 12912- 
12925 (2012), [http : //www . opticsinf obase ■ or g/oe/abstract . cf m?URI=oe-20- 12- 12912 

7. M. D. Feit and J. A. Fleck, "Computation of mode properties in optical fiber waveguides by a 
propagating beam method," Appl. Opt. 19, 1154-1164 (1980). 

8. M. D. Feit and J. A. Fleck, "Computation of mode eigenfunctions in graded-index optical fibers 
by the propagating beam method," Appl. Opt. 19, 2240-2246 (1980). 

9. G. P. Agrawal, Nonlinear fiber optics, second ed. (Academic Press, NY, 1995). 

10. D. Marcuse, Theory of dielectric optical waveguides, second ed. (Academic Press, NY, 1991). 

11. D. Marcuse, "Curvature loss formula for optical fibers," J. Opt. Soc. Am. 66, 216-220 (1976). 



12. K. D. Cole, "Steady-periodic Green's functions and thermal-measurement applications in rectan- 
gular coordinates," J. of Heat Trans. 128, 706-716 (2006); DOI: 10.1115/1.2194040 

13. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C, 
second ed. (Cambridge Univ. Press, NY, 1992). 

14. P. W. Milonni, The quantum vacuum: an introduction to quantum electronics (Academic Press, 
CA, 1994). 

15. M. Frigo, and S. G. Johnson, "FFTW Home Page," |http://www.fftw.org/| 



1. Introduction 

In earlier papers we described a physical processes that can lead to modal instability [1] 
in high power fiber amplifiers. We also presented simulation results that illustrated how 
periodic modulation of the pump or signal seed can dramatically reduce the instability 
threshold [2]. Although we described our model and methods in those papers we did 
not present the mathematical details. Here we do so. The individual parts of the model 
are all known, but the way they are combined to form a steady-periodic model of modal 
instability is original. 

First, we offer a brief reprise of the physics included in our fiber amplifier model. Mode 
instability is the degradation of output beam quality above a sharp power threshold [3] . 
Below threshold the light emerges in the fundamental mode, above threshold some is 
in higher order modes, usually LPn. Assuming most of the signal seed light is injected 
into the fundamental mode LPoi with a small amount populating a higher order mode, 
usually mode LPn, these two modes interfere along the length of the fiber. Because 
they have different propagation constants their interference creates a signal irradiance 
pattern that oscillates along the length of the fiber. Pump light is preferentially absorbed 
in regions of higher signal irradiance, and because a fraction of the absorbed pump is 
converted to heat, this creates a heating pattern that resembles the irradiance pattern. 
The heat pattern is converted to a similar temperature pattern, and the temperature 
pattern creates a refractive index pattern via the thermo optic effect. If the interference 
pattern is stationary, there is little or no phase shift between the thermally-induced 
index pattern and the irradiance pattern, so there is almost zero net transfer of power 
between modes. However, if the light in the higher order mode is slightly detuned in 
frequency from that in LPoi , the irradiance pattern moves along the fiber - down stream 
for a red detuning and up stream for a blue detuning. The temperature pattern also 
moves, but it lags the irradiance pattern, and this lag produces the phase shift necessary 
for power transfer between modes. Red detuning of LPn leads to power transfer from 
LPoi to LPn. The detuning that maximizes the mode coupling is set largely by the 
thermal diffusion time across the fiber core. This is approximately 1 ms, implying an 
optimum frequency detuning of approximately 1 kHz. 

When light in LPoi is transferred to LPn by deflection from the moving grating it 
experiences a frequency shift equal to the frequency offset, due to a Doppler effect pQ, 
so the transferred light adds coherently to the light already in LP 1 1 . This mechanism 
can cause a substantial transfer of the power from LPqi to LPn if the pump power 
exceeds a sharp mode instability threshold. This gain process can be categorized as 
near-forward stimulated thermal Rayleigh scattering. 

Our modeling approach is to develop as general a model as feasible, and our model 
includes all the physical effects just described. We use diffractivc beam propagation of 
a steady-periodic, time-dependent signal field. This field incorporates all fiber modes 
(including lossy modes) and all offset frequencies that are harmonics of the primary 
frequency offset. The time-dependent thermal profile is computed using either an 
alternating-dircction-implicit (ADI) integration method or a steady-periodic Green's 



function method. Mode coupling occurs through inclusion of the thermally-induced re- 
fractive index change in the index profile that is used to compute beam propagation. 
No analytic or semi-analytic expressions of mode coupling are needed, and coupling be- 
tween all modes is included. We use this approach because it permits the most accurate 
modeling, and because it makes adding various physical effects relatively straightfor- 
ward. The cost of such a general numerical model is long run-times, but using the 
methods described here one can model several meters of fiber per hour on a consumer 
grade desktop computer. 

Alternative mode coupling models based on similar physical effects have been pub- 
lished by other authors [3H5] but we will not review them here. 

2. The scalar, paraxial beam propagation equation 

Equation ([T]) is the scalar, paraxial beam propagation equation for an isotropic medium. 
It is derived in numerous papers on nonlinear fiber optics, for example [7|f5] . The variable 
E s is a complex envelope function that represents all spatial and temporal modulation of 
a monochromatic, plane carrier wave. The direction of propagation is z. This equation 
is appropriate for small index contrasts typical of step index, large mode area fiber 
amplifiers. 

dE s (x,y,z,t) i 2 i[k 2 (x,y,z,t)~k 2 ] 
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where k c = G) c n c iad/c is the wave number of the carrier wave in the cladding, and 
k(x,y,z,t) = 00 c n(x,y,z,t)/c. Here n(x,y,z,t) is the guiding index including the thermo- 
optic contribution, 

dn 

n(x,y,z,t) = n com (x,y,z,t) + —T(x,y,z,t), (2) 

where dn/dT is the thermo optic coefficient given in Table [T] The first term on the right 
hand side of Eq. ([TJ describes diffraction in a homogeneous medium with a refractive 
index equal to the cladding index n c . Here is the Laplacian operator in the transverse 
dimensions x and y. The second term adds a correction to account for the core guidance, 
including the usual core refractive index step plus the thermally induced index change. 
In general k(x,y,z,t) is imaginary to account for gain or loss, but we will use k(x,y,z,t) 
to represent only the real part, and write the much smaller imaginary part as a separate 
term with the coefficient g(x,y,z,t), 

dE(x,y,z,t) i j i[k 2 {x,y,z,t)-k 2 .] 

j; = 2^Vj_£(x,y,z,*) + - ^ L E(x,y,z,t)+g(x,y,z,t)E(x,y,z,t). (3) 

3. Integrating the propagation equation 

Equation ^ can be integrated by various methods - for example, using the split step 
method described below, or using finite-differences. Whichever method is used, the 
central issue is how to compute the thermally induced part of k{x,y,z,t) for use in the 
z integration. 

One approach is to integrate the field along the full length of the fiber at t = based 
on an initial temperature profile T(x,y,z,0). The temperature profile is then updated 
based on the heat deposited during the time increment At, plus the previous temperature 
profile and the thermal boundary conditions. The new temperature profile is used to 



integrate the field at time At along the full length of the fiber. This iteration of full z 
integration of E s alternating with stepping T by At is repeated for successive time steps 
to find the temperature and field histories, both over the domain (x,y,z,t). 

Alternatively, a time sequence of fields at z = can be used to compute T{x,y,t) at 
the input and that time varying temperature can be used to propagate the time varying 
field by Az. This is repeated for successive z steps. One limitation is that the temperature 
is solved in two dimensions rather than three so longitudinal heat flow is not accounted 
for. 

3.1. Steady-periodic condition 

Both of the methods just described can treat transient or steady-periodic cases. How- 
ever, our goal is to compute behavior only under steady-periodic conditions. We are not 
interested in starting from an initial transient and following the amplifier evolution to 
a steady state since that may require integration over a time longer than the thermal 
diffusion time from the core to the outer diameter of the fiber. The required settling 
time can be greater than 100 ms, perhaps much greater, depending on the fiber design. 

Because transient effects are of secondary interest, and also in order to achieve shorter 
run times, we base our model on the assumption that all transients have decayed and the 
only time dependence is periodic modulation of powers, mode content, temperature, etc. 
With this steady-periodic assumption we need integrate only over a single modulation 
period of approximately 1 ms. Wc first perform the time integration on the temperature 
equation to find T(x,y,t) for a single period at one z position, then we step the signal 
and pump fields by Az. 

The steady-periodic condition appears to correspond well with reported behavior of 
amplifiers operating near and slightly above their mode instability thresholds. Far above 
threshold the mode coupling may be more complicated [3]. The relevant time scale for 
the steady-periodic condition is the time for heat diffusion over the core rather than over 
the entire fiber. The temperature outside the core is determined by the time averaged 
heating in the core, but it cannot respond rapidly enough to follow the periodic heating 
variations within the core. 

3.2. Transverse heat flow approximation 

Our method does not allow longitudinal flow of heat. We justify this by noting the large 
difference in length scale for the transverse and longitudinal temperature variations. 
The important longitudinal length scale is the modal beat length. This varies with fiber 
design but it typically falls in the range 3-30 mm in large mode fibers. The transverse 
length scale is the core diameter which typically lies in the range 20-80 flm. The ratio 
of longitudinal to transverse length is of order 100:1. 

3.3. Narrow line width approximation 

Our model relics on a frequency offset between the light in different transverse modes. 
The frequency offset is usually in the range 300-3000 Hz. The assumption of periodic 
behavior implies the set of all allowed frequencies must be separated by multiples of the 
inverse modulation period. The light in any mode must consist of only these frequen- 
cies. The width of each frequency component is assumed to be much smaller than the 
frequency offset. This narrow line width is obviously not entirely realistic. Nevertheless, 
as we will discuss below in Sec. [TT1 this method can still be a highly accurate way of 
treating the problem of mode coupling in relatively narrow band amplifiers. 



3.4- Split-step integration in t and z 

Because it is only necessary to treat one oscillation period under the steady-periodic 
assumption, we perform the integration in z on a set of time samples of the field that are 
evenly spaced over one period. At the fiber input we specify E s (x,y,0,t) and the pump 
power. Using this information we compute the temperature profiles for each sampling 
time over the full period. We describe two methods for the integration of the heat 
equation below. The temperature profile for each time sample is then used to propagate 
the corresponding time sample of the signal field by Az. This is repeated until the end of 
fiber is reached. This method is used because it makes efficient use of limited computer 
resources, and because it can be made to run fast. 

We use a split-step, fast-Fourier-transform, beam propagation method (FFT BPM) 
to perform the z integrations. Split-step methods have been employed in numerical 
simulations of CW optical diffraction, including guiding in optical fibers, at least since 
the 1970s [TUB] . One period of the field is discrctized into N t time steps. We use the split- 
step method by applying it individually to each of the field time slices. Related FFT 
methods for dispersive propagation are also widely used to model propagation of short 
light pulses propagating in a single fiber mode [9]. Although we are also propagating 
pulses consisting of one beat time, it is best to think of the time slices as widely separated 
samples that are unaffected by dispersion. We will say more about this in Sec. [17] 

In the split-step method, Eq. ([3]) is used to advance each field time sample by Az in 
three steps. In the first step, linear diffractive propagation is applied to advance the 
field by Az/2 in a homogeneous medium with a refractive index equal to the cladding 
index. In this step, the beam propagation equation is used, keeping only the diffractive 
term on the right-hand-side 

dE s (x,y,z,t) i 2 

Tz = ^ v ^W,v). (4) 

In the second step, the field is advanced by Az without diffraction, keeping only the 
phases induced by the guiding and thermally induced index, along with the laser gain 
and linear loss contained in g(x,y,z,t). The propagation equation for this step becomes 

= i[eiXJ :f-^ E s ( X ,y,z,t) +8 ( X ,y,z,t)E s (x,y,z,t). (5) 
az 2k c 

The gain term will be described in more detail in the next section. In the third step, 
linear diffractive propagation is again applied to advance the field by Az/2. This method 
produces errors of order ^(Az 3 ). 

4. Algorithm 

These methods are incorporated in our algorithm which is diagrammed in Fig. [T] In the 
following sections, we describe each step of the algorithm in more detail. 

5. Read problem parameters and fiber parameters 

The first step in the algorithm is to read in the problem parameters. Table [1] lists 
them individually. Some of them are standard silica parameters. For example: thermal 
conductivity K — 1.38 W/m-K; specific heat C = 703 J/kg-K; density p = 2201 kg/m 3 ; 
thermo-optic coefficient dn/dT = 1.2 x 10~ 5 K _1 . Usually the pump wavelength is k p = 
976 nm, and the upper state ion lifetime is T = 901 /is. 
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Fig. 1. Flow diagram for split-step FFT propagation using alternating direction 
implicit integration (ADI) or steady-periodic Greens function method for tem- 
perature calculation. 



Table 1. Model input parameters 




We define our problem at each z- location on an (x,y,t) grid. The grid is equally spaced 
in (x,y,t), with number of points typically (N x = 64, N y = 64, Nt = 64). The spatial grid 
spans a domain of size L x x Ly with the core at its center, where L x and L y are typically 
«3x c/core- We choose as the thermal boundary condition a fixed temperature on the 
domain boundary. We assume there is a fixed frequency offset between modes equal to 
Acq. The period Y is defined as one cycle of the beat frequency, Y = 2n/Aco. The grid 
step sizes are then 



The parameters defining the doping profile NYb{x,y), the refractive index profile 
n C ore(x,y) and the linear loss cc(x,y) vary with (x,y) position. We usually use super- 
Gaussian profiles for these, with super- Gaussian coefficient of order 40. The FWHM 
values for the super-Gaussian diameters are d C oie and c/dopant- The pump is confined to 
the pump cladding of diameter c/ c iad- The absorption and emission cross sections for 
pump and signal are G°, O e p , a", o~f . 

The forward and backward time-averaged input pump powers, P^ and Pg, may be 
periodically modulated by specifying a pump modulation function M p (t). Similarly, each 
signal mode LP mi „ with time averaged power Pf nn can be modulated by specifying its 
modulation function M s m „{t). The J mn are mode specific losses (non heating). 

The fiber length is L, and the fiber is bent to a radius of i^bend- The integration step size 
Az is typically set to a few microns. We store information about the field every A samp i e 
integration steps. Typically the stored information includes the local time sequences for 
the mode content, total signal power, effective area of the signal, and pump power. 

6. Calculate modes 

The next step of the algorithm is to compute the signal mode profiles. Here, we present 
a general method capable of handling bent fiber and arbitrary refractive index profiles 
of low contrast. This procedure is from Marcuse [TO]. It is not a vector method so it is 
not appropriate for air holes or other guiding structure with high index contrast. Vector 
adaptations that permit high index contrast are available, but we do not discuss them 
here. 
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On the rectangular grid (0 < x < L x ) by (0 < y < L y ), the functions 



with n and v being integers, form a complete, orthonormal set, normalized so their 
areal integral is unity. Any transverse mode \j/(x,y,z) can be written 

y(x,y,z) = e-^ £ £ C MV ^ v (x,y) (10) 
li=l v=l 

where j3 is the propagation constant for that mode, C^ v is the mode's eigenvector 
(coefficients of the basis set 0^ v from Eq. ©) and \j/ is a solution of the Hclmholtz 
equation 

d 2 w d 2 w d 2 w , , , , 

-^r + -^ + -^r+ n 2 (x,y)%v = o (n) 

where £ = ffl c /c. We set an upper limit on the number of sin-sin terms to include in 
the expansion of y/ (M* and My), and substituting this expansion into the Helmholtz 
equation yields 



M t My 

/x=l v=l 



2 2 



fyv(x,y)=o. (12) 



Next multiply this equation by and integrate it over x and y. We define n e ff = /k 
and we divide the equation by k 2 and add and subtract n^ lad under the summation sign 
to obtain 

m x My 

£ £^Vv', jiv C MV = («eff- w clad)C M 'v'- (13) 
|i = l v=l 



The A matrix is 



ckdyj [« 2 (x,y) -n^iad] ( Vv'(*)) ; )0aiv( x >) ; )} 



Ly rLy 

io 

-f(| + Z|) 5 »' 5 -" (14) 

Equation (|13p is an eigenvalue equation of the form 

AC = n 2 C (15) 

to be solved for eigenvalues h and eigenvectors C. The solution eigenvectors are substi- 
tuted in Eq. (|10l) to compute the cigenmodes. Then the effective refractive index « e g- of 
mode LP„ V , is related to its eigenvalue n mj „ by 



n eS (m,n) = \jn 2 n +n\ A& . (16) 

The beat length between LPoi and LPn can be found using their values for rc e ff. 

For simple guiding index shapes such as a step index in unbent fiber, the integrals 
in Eq. (|14[) can be integrated analytically, otherwise they can be evaluated numerically. 
We use numerical integration because of its generality. 

The modes computed this way are the low power modes that do not take into account 
a nonuniform temperature or any other irradiance-dependent index changes. These low 
power basis modes will be used to construct input fields and to analyze propagating 
fields into modes. 



6.1. Refractive index profile 

The method described is general enough to consider arbitrarily-shaped refractive index 
profiles. For simplicity, we typically use a top hat-like two-dimensional super- Gaussian 
of high order to construct our refractive index and doping profiles. The super-Gaussian 
index profile is computed using 



n(x,y) 



: ^clad H~ (^core ^clad) 

xexp ^-ln2|[2(x-x )/c/ C ore] 2 + [2{y - yo)/d coie ] 2 } ^ 



(17) 

where n c i a d and 

^core the refractive indices of the cladding and core, (jcojJo) are the 
coordinates of the core center, <i core is the core diameter, and S is the super- Gaussian 
coefficient (we typically use S = 40). 

We approximate the effects of bending the fiber by adding a refractive index ramp 
Abend to n(x,y) in the plane of the bend so the index is higher on the outside of the bend. 
If bending is in the x plane, this added ramp can be written 

"bendC*^) =n(x,y)(x-x )/R hend . (18) 

Bend losses are calculated using the method of Marcuse [IT] . In the numerical model 
bend and other losses from the core are enforced by included an absorbing layer near 
the grid boundary. 

6.2. Normalization 

In order to easily relate our calculated modes to power in watts, we introduce a nor- 
malization factor Nm y n for each mode constructed from Eq. (|10p , where N m>n is the value 
which satisfies 

££o r 
2 

The normalized mode M m „(x,y) is defined by 

u m ,n{x,y) =N m>n \lf m>n (x,y). (20) 



J dxdy n(x,y) N m , n y m , n (x,y) =IW (19) 



7. Set up temperature solver 

Next the temperature profile must be solved for each time point in the period Y using 
the periodic heat source Q. In an isotropic homogeneous medium, the heat equation is 

where T = T(x,y,t) is the temperature, Q — Q(x,y,t) is the heating source, K is the scalar 
thermal conductivity, p is the density and C is the specific heat capacity. No d 2 /dz 2 
term is included in Eq. (f2"Tj) because we don't allow longitudinal heat flow. We consider 
only the thermal boundary condition of a fixed temperature on the (x,y) grid boundary. 
Alternative boundary conditions can be dealt with by modifying our procedure, but we 
do not discuss them here. 

We will describe two methods of solving Eq. (|2"T|) for T. The first method, the Green's 
function method for steady-periodic heating, involves computing the temperature over 
the (x,y) grid as a sum of independent contributions, one from each heated grid point. 
The second method, the alternating direction implicit (ADI) integration, involves step- 
ping in time via successive matrix multiplications involving all pixels contributing to- 
gether. 



7.1. Calculate Green's functions for steady-periodic heating 

To use the Green's function method we must first compute the steady-periodic Green's 
function which gives the temperature contributions over the entire (x,y) grid due to 
the periodically heated point (x 1 ,y'). These functions describe the temperature rise due 
to a steady-periodic heat source at a single modulation frequency. Formulations of the 
Green's functions for different boundary conditions arc detailed in |12j . For temperatures 
clamped at all four sides of the grid, the complex valued Green's function is 

G(x,y,ft)|x',/) = f.F n {y,y l )P n {x,x',CO) (22) 

n=0 

F„(v,/) = ^Lsin^v) S in(fy) (23) 

P„{x,x',a>) = | exp [ — a„(2H- \x-x'\)] — evp[—O n (2H—x—x')] 

+ exp [ - a„ \x -x'\] - exp [ - o„ (x + x') 
a n (l -exp[-2tr„H]j 

where 

ft) is the frequency of the heat, (0 <x < H), and (0 < y < W). 

We compute these functions for several harmonics, making sure the time grid is 
sufficient to resolve them. That is, we find G(x,y,co,„\x' ,y') for ft),,, = (0,1,2,3,...) x Aft) 
where Aft) is the specified frequency offset. The use of this set of Green's functions to 
compute T(x,y,t) is described in Sec. 113. T1 

7.2. Calculate ADI matrices 

An alternative, more general method of solving the heat equation is the alternating 
direction implicit (ADI) method |13j . This method does not require periodicity. When 
it is applied to a steady-periodic problem, it requires many iterations to match the 
steady-periodic requirement, with each iteration operating on the previous iteration's 
data. This sequential process means a simple parallel computing scheme is difficult to 
implement for ADI integration. 

The derivation of the equations for this method begins by using finite difference 
definitions of the partial derivatives in Eq. (|21[) . The numerical integration of the heat 
equation by time step At is split into two halves. We first integrate explicitly in y and 
implicitly in x for one half time step At/2, then integrate explicitly in x and implicitly 
in y for another half time step At/2. 

We start with the combined expression of implicit-x/explicit-y integration 

T r+Af/2 _ T t i T t+te/2 _ 1T t+At/2 T r+Af/2\ %oc 

lx ,y l x,y — Vx+Ax,y llx -y + l x-Ax,y ) ~T 



Ay 

~2 

At 



Tx,y+Ay _ 2T* >y + T*y_ Ay ) — (25) 
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or, rewritten in matrix form, 
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x ' 



(26) 
(27) 



where A x is the matrix for implicit integration in the x-direction, and B y is the matrix 
for explicit integration in the y-direction. These matrices are tridiagonal, and A x has 
size (N x — 2) x (N v — 2) while B y has size (N y - 2) x (N y — 2). The matrices are 



A x = 
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(28) 



(29) 



Matrix A, has diagonal elements (1 +At) and off-diagonal elements —At/2, while matrix 
B y has diagonal elements (1 — A,) and off-diagonal elements Ay/2, where 



KAt 
KAt 



J* = -^72 (3°) 



h = ^72- (31) 



In the second half time step, implicit/explicit integration order is reversed 

T t+At _ T t+Ai/2 _ / T r+Af/2_ 9T r+Af/2 T r+Ar/2\A r 
l x,y 1 x,y — \ l x+Ax,y ll x.y + 1 x -Ax,y ) 2 

.(ft+At _- >T t+At , T i+At fnn\ 
~T~Vx,y+Ay Zi jy> + I x.y-Ay) 2 ^ > 

Af gf+3A</4 



2pC 



Or, rewritten in matrix form, 

AT t+At = p+At/2 B + _^_Qt+3At/4 

jt+M = 4 -l r f+Af/2 2?je+ ^L^l fi f+3Af/4 j (33) 

where A y and B x follow a definition similar to Eqs. 1281 and 1291 with X x swapped with A, 
and redimensioned as appropriate. 



Combining the two half time step integrations, using the ADI to integrate by one full 
time step Af becomes 

r' +Af =AJ%T t A- l B x +^-A-' (V +A ' /4 A; 1 fi.v + e r+3Ar/4 Y (34) 

where, because we only have Q defined on integral numbers of time steps, we use linear 
interpolation to find Q t+At / 4 anc j g'+3A;/4 f rom t ne heat at t and t +Af. This method is 
converged to (3 {At 2 ). 

8. Construct propagation phase array 

Linear diffractive propagation that comprises the first and third steps in the split-step 
propagation is best evaluated in (k x ,k y ,z,t) space. Equation ^ can be transformed to 
this space by inserting the following form of E s (x,y,z,t) 

E s (x,y,z,t) = ^ J J E s (k x ,k y ,z,t)e ik " x e ik yydk x dk y . (35) 

The transformed equation is 



dE s (k x ,ky,Z,t) . kl - A 

^ = -i^E s (k x ,ky,z,t)-i^E s (k x ,ky,z,t). (36) 



Advancing the field E(k x ,k y ,z,t) by Az/2 consists of shifting the phase of each (k x ,k y ) 
plane wave component by 

t(k x ,ky) = -^ y -^-^. (37) 

The inverse two-dimensional Fourier transform is used to convert E s (k x ,k y ,z,t) back to 
E s (x,y,z,t ). In practice, fast Fourier transforms (FFTs) are used to efficiently convert 
fields between (jt,y,z,f) and (k Xl k y ,z,t) spaces. 

9. Construct input signal field and pump powers 

9.1. Signal 

We use the set of normalized modes M mj „(x,y) defined in Eq. ([20]) to construct the input 
signal field 

E s (x,y,z,t)\ = £ Jpj~ n u m , n {x,y)M^ n {t) (38) 

modes 

where P£ n is the power in the (m,ra) th mode and M* ln (t) the periodic modulation 
function for that mode. For example, to model an amplifier with LPoi populated 
by unshiftcd light, and LPn populated by light detuned by A©, we set Mq l (t) = 1, 
M\ j(f) = exp(— zAfiM). If instead we wish to amplitude modulate the signal in both 
modes with a simple sinusoidal envelope, we make each modulation function M'l in {t) = 
1 +0.25acos(Aof) where a is the small peak-to-peak power modulation. We can also 
use more complicated modulations, as long as they are periodic. For more complicated 
modulation we must include an adequate number of harmonics in the Green's function 
temperature solver. 



9.2. Pump 

Defining a co-propagating pump input is simple, but since we start the integration at 
the signal input end, specifying the counter-propagating pump that gives the target 
pump at the opposite end can be more difficult, especially if amplitude modulations 
are involved. To keep the two pumps separate we define two quantities P?(z,t)\ z= o and 
P£(z,t)\ z= o f or forward- and backward-propagating pumps. 

To generate modulated input pump powers, we follow a technique similar to the 
signal modulation considered above for each pump. The model is capable of treating 
an arbitrary periodic pump modulation, provided we include an adequate number of 
harmonics if solving the heat equation using the Green's function method. 

10. Propagate signal field 

To reiterate the propagation described earlier, the signal is propagated by Fourier trans- 
forming each time slice of the field E s (x,y,z,t) to E s (k x ,k y ,z,t) using 2D-FFTs. The phase 
of each plane wave component is advanced by the propagation phase <j)(k x ,k y ) given in 
Eq. (f37j) . and the inverse 2D-FFTs of E s (k x ,k y ,z,t) gives the updated field E s (x,y,z,t). If 
the propagation step is not the first or last half-step in the fiber, the two consecutive 
half-step propagations are combined by advancing the phase by 2<j)(k x ,k y ). 

11. Update the field for laser gain 

In the non diffractive part of the split-step procedure we update the signal field by 
adding the guiding and thermally induced phases plus any laser gain and linear loss. To 
compute the gain we find the upper state population profile in (x,y,t) from the signal 
irradiance profile, pump power, and doping profile, and use it to compute the signal 
gain and pump loss. The (steady-state) upper state population is 

, A _ P p o? } /h Vl A P + W/hv s 

nu[X,y,t} P p {o* + o<)/hv p A p +I,(o' + o>)/hv, + l/T 1 ' 

where I s = I s (x,y,t) is the signal irradiance, V s and V p are the signal and pump frequen- 
cies, of and o/ are the signal absorption and emission cross sections, C7 ; " and G e p are 
the pump absorption and emission cross sections, T is the ion upper-state lifetime, P p 
is the pump power, and A p is the area of the pump cladding. We assume that the cross 
sections and lifetime are independent of temperature. The effective ion lifetime at typi- 
cal amplifier powers is a few micro seconds so at modulation frequencies of order 1 kHz 
the steady state solution is appropriate. 
The change in signal field is given by 



dE s (x,y,t) _ 1 
dz 2 



Nn(x t y)E s (x,y,t) - ^a(x,y)E s (x,y,t), (40) 



where Nyb{x,y) is the Yb 3+ doping profile, a(x,y) is the linear absorption coefficient. 
This method is general enough to treat an arbitrarily-shaped doping profile, but we 
typically consider super-Gaussian doping profiles similar to the one defined in Eq. (|17|) . 
The linear absorption coefficient Ct can be non- uniform in (x,y) to accommodate a 
photodarkening model. We assume that all the power absorbed due to <x(x,y) is turned 
into heat. Referring to Eq. (JSJ) , the laser gain and loss term g(x,y,t)E s (x,y,t) is given by 
the right hand side of Eq. (|i0"|). 



12. Update pump powers 



We include both forward- and backward-going pumps, with the total pump power given 
by 

(41) 



P P =P P f 



P p is assumed to be uniformly distributed across the pump cladding. The change in the 
pump power is computed directly from the ion inversion, rather than from the signal 
increment. This allows us to include linear signal loss and fluorescence loss correctly. 
The total change in pump power is given by 



dPp 
dz 



[o a p + e p )n u {x,y) - 0° NYb(x,y)dxdy. 



The loss is apportioned between the forward and backward pumps according to 

dP/ &P p P p f 



dz 



dz Pn 



(42) 

(43) 
(44) 



dz dz P p 
13. Compute T(x,y,t) and thermo-optic index 

For the computation of T(x,y,t) we use either the Green's function method or the 
ADI method. We calculate the heat deposition rate from the absorbed pump and the 
quantum defect according to 



Q(x,y,t) =JV Yb (x,y) 



oZ-{oi + o e p )n„(x,y,t) 



Ppif) 



+ a(x,y)I s (x,y,t), (45) 



where the upper state population n u (x,y,t) is given by Eq. (|39l) . 
13.1. Green's function method 

The heat as a function of time over the full cycle at each (x',y') pixel is resolved into its 
Fourier components CO m = (0, 1,2,3...) x Act) by performing a temporal Fourier transform 
on Q(x',y',t) 

N,-l 

q(x J ,y',(O m )=AxAy £ Q(x',y',t i ) exp (-;©,„?,). (46) 

1=0 

Here, q(x' ,y' ,(O m ) is a complex coefficient that includes the phase as well as the ampli- 
tude of the heat deposition. These q(x',y',(O m ) values are used to weight the Green's 
functions in computing the steady-periodic temperature over the entire transverse grid. 
We always include frequencies (0, 1) x Act), and we add higher frequency terms as needed 
for convergence. The time grid must be fine enough to resolve the highest frequency. 
Higher frequency terms are usually necessary only above the mode instability threshold. 

Using the Green's functions computed as described in Sec. 17. 1[ the temperature 
T(x,y,t) is found by summing the contributions of each q(x' ,y' ,co m ) according to 



T(x,y,t) = Real 



X^O^'y^m) G(x,y,(O m \x',y') exp(i(0,„t) 



(47) 



13.2. ADI Method 



In Sec. 17.21 we described a method of integrating the heat equation using the ADI 
method. This method uses Eq. to integrate one time step of At. For steady-periodic 
heating, we enforce the steady-state criterion by integrating for several periods, reusing 
the heat Q(x,y,t) from one period to the next. We terminate this process when the 
difference in temperatures between periods has reached an acceptably small residual. 

The ADI method is more general than the steady-periodic Green's function method, 
because it does not require periodicity. It can be used to model transient heating, for 
example, if the steady-periodic condition is not enforced. It is unconditionally stable, 
but its accuracy suffers if the time-step is too large. 

14. Add thermo-optic and guiding index phases to the field 

In the non diffraction portion of the split-step propagation we also advance the phase 
of the field to account for the guiding index plus the thermal index over Az using 

Hx,y,t)=Az ^ (x ^l ) ~ k2c . (48) 

The phase can be rewritten using Eq. @ as 

At a a ac i \ t a i , H*,;y,0-"ciad] 2 N \ 

(j>(x,y,t) = Az— \ [n(x,y,t)-n daii H . (49) 

c y 2« clad J 

We write it in this form merely to show that the phase is (co c AzAn/c) plus a small 
correction from the quadratic term. 

15. Resolve time-dependent modal power contents 

We use the normalized modes defined in Eq. (|20|) to resolve the content of the signal 
field E s (x,y,z,t) into fiber modes. We compute the inner product of the field and the 
normalized mode from 

J7 I A ffdxdy E s (x,y,z,t)u m .„(x,y) 

ffl '" (z ' } = — — r~\T 2 — ' (50) 

JJ dtdy [u m ,„[x,y)\ 

where F nh „(z,t) is a complex quantity. The mode amplitude F mfl (z,t) is converted to 
power in watts using 

2 

Pm,n{z,t) = F m>n ( Z ,t) . (51) 

15.1. Spectral analysis 

To find the frequency spectrum of mode (m,n) we Fourier transform F mM (z,t) from time 
to frequency. The allowed frequencies for the periodic function F m ^ n (z,t) are the co m . 

15.2. Compute effective area 

We also compute the effective area of the signal field using 

2 



A eS (z,t) 



Jf\E s (x,y,z,t)\ 2 dxdy 



ff\E s (x,y,z,t)\ 4 dxdy 



(52) 



16. An example computation 

In Fig. [2] we show a sample result produced by the model. The surface shows the power 
in LP 11 over one LPoi-LPn beat length at the input end of an amplifier versus t and 
z. The motion of the surface ridges reflect the change in phase of LPn relative to LPoi 
due to the difference in propagation constants for the two modes. In this example the 
LPn seed light is Stokes shifted by 600 Hz relative to LPoi, but the details of the 
fiber are unimportant for this illustration since all large mode fiber amplifiers produce 
qualitatively similar surfaces. A small fraction of the gain is due to laser gain but most 
is due to thermally induced mode coupling. 

The offset frequency that produces the highest mode coupling gain varies along the 
length of the fiber due to the changing shape of the heat profile. We integrate a set 
of frequencies and powers over the full length of the fiber to find the lowest threshold. 
When operating near or below threshold the gain curve has a well defined frequency of 
highest gain even though the shape of the gain varies somewhat along the fiber due to 
changing heat profiles. 




Fig. 2. Power in LPn versus time (T) and distance (Z) at the input end of the 
amplifier similar to the one specified in [2]. This illustrates the growth of the 
low power mode over one beat cycle and over one beat length. The pump is 
unmodulated in this example. 



17. A narrow band width model for broad band light? 

Our model assumes that inter modal frequency shifts of order 1 kHz are responsible 
for thermal mode coupling. How is it possible that a model based on such small shifts 
can be applicable to light that has much broader line widths? Few seed lasers have sub 
kHz line widths, and in fact, the seed light is often intentionally phase modulated to 
widths of several GHz to combat SBS. Our answer is that even broad band light can 
be frequency shifted as a whole by 1 kHz. Further, for phase modulated light the beat 
between two waves with identical, but frequency shifted spectra is just as strong as the 
beat between two frequency shifted monochromatic waves. Additionally, the frequency 
shift in our model is caused by coupling two modes via a moving thermal grating. 
The Doppler frequency shift induced by this process shifts the entire spectrum of the 
scattered light so the interference between modes has full strength. This means that 
both in the mode coupling process and in the mode beating process responsible for 
heating, light with a broad spectrum behaves the same as monochromatic light. Also, 
a phase modulation of several GHz has little impact on the ion population evolution. 
The steady state expression for the upper state population is still accurate. 

This equivalence of narrow and broad band light breaks down if the two interfering 
fields shift out of time synchronization. This can happen for large band widths because 
different modes have slightly different group velocities. The difference in modal propa- 
gation times is of order 1 ps/m, so an estimate of the line width limit is 100 GHz for a 
typical fiber amplifier. However, this is substantially larger than typical SBS suppressing 
line widths. 

If the signal light consists of a train of short pulses, and the pulse separation is shorter 
than a few micro seconds, the ion population responds primarily to the average power. 
Slow modulations of the average power, in the kHz range, arc still effective in driving 
the thermal grating that leads to mode coupling. If the pulse width is more than a ps 
or so the pulses in different modes will stay synchronized well enough to produce strong 
inter modal interference. 

18. Sources of frequency offset HOM seed light 

It is likely that in most amplifiers an amplitude or spectral modulation of the pump 
light [2], combined with accidental injection of a small fraction of seed light into LPn, 
produces the initial frequency shifted light in LPn. Alternatively, amplitude modu- 
lation of the seed, plus accidental injection of a fraction into LPn could provide the 
initial light. If both of these modulations can be sufficiently suppressed, the initial light 
would be provided by quantum noise. This unavoidable initial light would produce the 
highest possible instability threshold. Our model does not incorporate a true quantum 
noise model. Instead we estimate the starting noise power from a classical stochastic 
electrodynamics noise model |14j often used in estimating thresholds for lasing or for 
Raman gain. The noise level is set to half a photon per mode. In our case there is a 
single transverse mode and a gain bandwidth of approximately 1 kHz. The expression 
for starting power is the same as for Raman noise power, 

P = hvAv = (6.6 x 10~ 34 )(2.8 x 10 14 )Av = 1.8 x 10~ 19 Av (53) 

which gives 10~ 16 W for a bandwidth of 0.5 kHz and a wavelength of 1060 nm. This 
is the minimum starting power in the higher order mode. The gain process will pick 
out from the sea of quantum noise the light that best matches the light in LPoi both 
in amplitude and phase over the beat cycle. This is the light with highest gain. There 



will be fluctuations in frequency and power, but the average power within the gain 
bandwidth will be well approximated by Eq. (j53jl . 

19. Desktop computer implementation 

19.1. Memory requirements 

Storing the full, four-dimensional double-precision, complex [64 x 64 x 64 x L/Az] array 
of the signal field would require on the order of 1 — 10 terabytes in a meter long amplifier 
with step sizes of a few microns. Therefore, we do not store E fields at each step. We only 
store computed properties such as modal content F m ^,(z,t), total signal power I(z,t), and 
effective area A e ff(z,f) at positions separated by (A samp i e x Az). This reduces the amount 
of memory required to at most a few gigabytes. 

19.2. Parallel computing 

A substantial portion of the model's run time is spent solving the thermal problem. 
This makes the Green's function method attractive because it is generally much faster 
than the ADI method. One reason the ADI is relatively slow is that it is necessary 
to integrate through several cycles of the heat to ensure steady-periodic condition is 
enforced. 

Equally important, the Green's function method is easy to parallelize while the ADI 
method is not. The ADI method is sequential by nature. It requires updating the entire 
grid to advance the time by At. In contrast, the Green's function involves a summation 
of the temperature contribution from each heated pixel, so the calculations for pixels 
are independent and easy to run in parallel. This parallelization is relatively simple to 
implement using the shared-memory multiprocessing library OpenMP. 

In addition, we can parallelize other steps of the model. Propagation of the signal 
field array can be performed independently for each time slice of the field. Similarly, 
the laser gain calculation, modal content decomposition, and effective area calculation 
can be performed independently for each of the time slices. 

19.3. Execution speed considerations 

We have both MATLAB and Fortran versions of our model. Both versions use the FFT 
library FFTW [T5]. Our experience is that FFTW is substantially faster than other 
FFT routines in Fortran. 

OpenMP, the shared-memory multiprocessing library and compiler directives, is im- 
plemented in the version of Fortran we use, GNU Fortran 4.6.3. Using a desktop com- 
puter based on an Intel Core-i7 3770 (Ivy-Bridge) processor with four physical cores, 
we are able to obtain an excellent speed up. Using four threads instead of one, the time 
required to run the same fiber setup decreases by a factor of « 3.2. 

Another important speedup was obtained by solving the thermal problem and laser 
gain equations on a coarser z-grid than the FFT propagation problem. This can dra- 
matically reduce the model run time as well. Care in choosing this grid is important 
because a grid that is too coarse can reduce the computed gain and increase the insta- 
bility threshold power. 

The run times on our computer lie in the range of 0.25-1.5 hour/m depending mostly 
on the z step size and the number of harmonics included in the Green's function. Larger 
cores use larger z steps, and near-threshold runs require only a single harmonic, per- 
mitting times near 0.25 hour/m. 



20. Approximations of model 

All numerical models make judicious approximations. A brief list of ours follows: 
Single X pump with single absorption and emission cross sections 
Pump power is uniformly distributed across the pump cladding 
All signal light is identically polarized 

Steady-periodic heating is required for application of Green's function 

Fixed period Y, which allows only signal frequency offsets 1/Y x (0,±1,±2, ...) 

Thermal boundary condition is fixed temperature on square boundary 

Thermal boundary size approximately three times core diameter 

Heat equation solved in two dimensions (no longitudinal heat flow) 

Thermal properties assumed uniform and isotropic 

Temperature dependence of cross-sections, dn/dT, and T not included 

No refractive index dependence on n u 

Signal bandwidth < 100 GHz 

Low contrast refractive index profile, e.g. no air-holes as in PCF 

Upper state population n u follows I s instantaneously (steady-state expression) 

21. Attributes of model 

Attributes of our model include: 

Highly numeric - general and simple to add additional physical effects 

Model a variety of refractive index, linear absorption, and doping profiles 

Steady-periodic eliminates long integration times before steady state 

Steady-periodic model produces well defined thresholds 

The ADI method can be used to study transient behavior if desired 

All transverse modes are automatically included 

Thermal lensing automatically included 

Comparatively short run times and minimal memory requirements 
Variety of thermal boundary conditions possible using Green's functions 
Green's function method offers large speedup using multiple processors 
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